Method to process images obtained by tomography or few view tomosynthesis

ABSTRACT

The invention concerns a method for processing images obtained by tomography or tomosynthesis, the method comprising: acquiring a plurality of 2D projection images of an object using an imaging system, the acquiring being defined by Rf=s, wherein s is a vector of the acquired projections, R is a projection operator which models the imaging system and f is the 3D image of the object to be reconstructed with knowledge of R and s; and processing the acquired 2D projection images, wherein the processing comprises applying an iterative algorithm defined by its iteration, wherein at least one set of processed projection images is generated at each iteration, the iteration being defined so that at each iteration the 3D image of the object is a linear function of the processed 2D projection images according to the property f (n) =  R p (n) , wherein  R  is a matrix such that  R Rf is a 3D image.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The invention concerns the reconstruction of a two-dimensional (2D) orthree-dimensional (3D) image of an object, for example part of a regionof interest in a patient, on the basis of a set of one-dimensional ortwo-dimensional views respectively of the region of interest, taken fromdifferent positions by an imaging system around the region of interest.The invention finds particular application in medical imaging bytomography reconstruction or Few-View Tomography.

2. Description of the Prior Art

With tomography it is possible to produce slice images of a region ofinterest of an object, by acquiring projections.

FIG. 1 schematically illustrates the acquisition of 2D images of a bodyorgan and the reconstruction of this organ as a 3D image bytomosynthesis.

X-rays R derived from a source S are emitted at different angles (1, . .. , i, . . . , n) towards the organ O. After passing through the organthey are detected by the detector D_(ct) forming a set of 2Dprojections, s₁, . . . , s_(i), . . . , s_(n). It is to be noted thatthere are as many acquired 2D projections as there are angles underconsideration. If all the angles cover at least a semicircle, the termtomography is used. When this condition cannot be met (for example: arcof a circle of less than 180 degrees, straight line) the termtomosynthesis is used.

Acquisition is obtained using a detector positioned facing an X-raysource, e.g. a digital camera.

One application of tomography is the detection and characterization of alesion in an organ, e.g. a cancer lesion. For the breast, tomosynthesisis used.

The acquired 2D projections are used to reconstruct a 3D image of theobject. This 3D image is more precisely a 3D map of the X-rayattenuation coefficients of the medium through which the rays havepassed.

It is this mapping that is used by a radiologist to interpret this imagein relation to the observed differences in contrast.

3D reconstruction of the image is a costly, complex operation, inparticular if it uses an iterative algorithm.

An iterative 3D reconstruction method is known. This method is based ona discrete, matrix expression of the topographic reconstruction problem.

More precisely the problem can be modelled by the following equation:

Rf=s

in which s is a vector of the acquired projections, R is a projectionoperator which models the topographic imaging system and f is the 3Dimage of the object to be reconstructed.

The problem to be solved with tomographic reconstruction is to determinef having knowledge of s and R.

Exact reconstruction of the 3D image f in theory requires a set ofmeasurements at least as large as the number of unknowns forming the 3Dimage, measured under the constraints of tomography i.e. measurementswith a very small sampling pitch and a rotation of 180 degrees about theobject.

In practice, this may prove to be impossible for the following reasons:

Firstly, in examinations of tomosynthesis type, it is not possible toobtain measurements at 180° around the object to be imaged, and secondlythe smaller the sampling pitch the higher the number of measurements,which increases the duration of examinations and hence the applied X-raydose.

Therefore the case is frequently encountered in tomography in which thesize of the 3 image is much larger than the size of the data, whichimplies solving a system having a very high number of possiblesolutions.

The choice of reconstruction method allows a 3D image to be found thatis an acceptable solution for the intended application.

Additionally, known space reconstruction methods of 3 images are in theform of an iterative algorithm defined by its iteration A whichgenerates a reconstruction f^((n)) of the solution from a reconstructionf^((n-1)) such that:

f^((n))=A[f^((n-1))].

Starting with an initial reconstruction denoted f⁽⁰⁾ in the series ofreconstructions f^((n)), the reconstruction at iteration N is obtainedas:

f^((N)) = A[f^((N − 1))] = A[A[f^((N − 2))]] = A²[f^((N − 2))] = … = A^(i)[f^((N − i))] = … = A^(N)[f⁽⁰⁾].

This algorithm requires the storing of initialization steps, ofintermediate reconstructions and archiving of the reconstruction, i.e.one or more 3D images depending on whether storing is madesimultaneously or successively, and represents a generated amount ofdata that is much greater than the amount of acquired data (low numberof projections).

Some known methods use discrete, matrix (hence linear) operations. Someare associated with solving the linear system Rf=s, others areassociated with minimization of the quadratic error ∥Rf−s∥² sir and ofthe associated system R′Rf=R′s in which R′ is the transposed matrix ofR, others generalize the preceding approach by solving a system RRf= Rsin which R is a matrix such that RRf is a 3D image. Reference can bemade to the document <<CIARLET, P. G., Introduction á l'AnalyseNumerique MairieieIle el á l'Optinnsation, Masson. Paris, 1982>>.

More recently, it was proposed to calculate a filter iteratively, andonce and for all, the filter being obtained by iterative reconstructionof a test object and to apply this filter to every object to be imagedat the acquired projections in order to obtain the reconstructed volumeof the object without any iterative computing. In this respect,reference can be made to document U.S. Pat. No. 7,447,295. In this way,the required computing power is reduced compared with iterative methods.However, the substitution of iterative reconstruction by filteringgenerates loss of data. It is therefore a near-approach method.

SUMMARY OF THE INVENTION

The invention proposes an iterative method which does not necessitatethe storage of intermediate 3D images.

According to a first aspect, the invention concerns a method to processimages obtained by tomography or tomosynthesis comprising: anacquisition step to acquire a plurality of 2D projection images of anobject, acquisition being defined by Rf=s in which s is a vector of theacquired projections, R is a projection operator which models theimaging system and f is the 3D image of the object to be reconstructedhaving knowledge of R and s: a processing step to process the acquired2D projection images.

The method according to the first aspect of the invention ischaracterized in that processing of the projection images consists ofapplying an iterative algorithm defined by its iteration V which at eachiteration, n=1, . . . , N generates at least one set of processedprojection images p^((n)).

Iteration V is defined so that at each iteration the 3D image of theobject is a linear function of the processed 2D projection imagesaccording to the property f^((n))= Rp^((n)), in which R is a matrix suchthat RRf is a 3D image.

Having regard to the constraint associated with the iterative algorithmV, only one set of 2D projection images is handled. This avoids, at anytime in the method (initialization n=0, intermediate computing 0<n<N,archiving n=N) having to store the entirety of a 3D image whose size ismuch larger than the set of acquired measurements.

Also, the equation f^((N)))= Rp^((N)) is determined in sufficientlysimple manner so as only to be calculated during the viewing step. Inparticular, this equation requires all the processed 2D images p^((N)),but only the computing and storing of that part of f^((N)) which isviewed.

Since it is neither possible nor necessary to view all the sectionplanes of a volume simultaneously, but only successively, it is nevernecessary to use storage means covering the entire size of the 3D imagef^((N)). To summarize, the storing of an entire 3D image is necessaryneither during the reconstruction step, nor during archiving, nor duringthe viewing step.

For known algorithms defined by their iteration A associated with thematrix RR in the 3D image domain which therefore requires the handling,storing and archiving of at least one stored 3D image in its entirety,it is shown that there exists an algorithm defined by iteration V whichgenerates a sequence of images processed in the form of 2D projectionimages and a linear constraint R allowing the generation of areconstruction of the object from processed projections, such that thereconstructions obtained by V and R are mathematically equivalent tothose generated by A, without it ever being necessary to store theentirety of a 3D image.

Also, the computing time is proportional to the sampling of each 3Dimage to be reconstructed, whereas the size occupied by storing 2Dprojection images is constant.

It is therefore possible to adjust the sampling of the reconstructed 3Dimage beyond the storage capacities of a processing unit underconsideration, and therefore to implement on this unit an iteration Vequivalent to an iteration A that cannot be implemented on this unit butonly on another, considerably increased storage unit.

Other characteristics of the method according to the first aspect of theinvention are the following:

the method comprises a 3D reconstruction step (S2) of the object to bedetermined f^((N))= Rp^((N)) in which N is the number of fixediterations;

at each iteration n=1, . . . , N, the algorithm V generates a set ofprocessed projection images p^((n))=V[p^((n-1))], iteration V beingdefined for any vector p in space of the 2D projection images in thefollowing manner V[p]=(1−ρR R)p+ρs in which p⁽⁰⁾=ρs or else p⁽⁰⁾=0, s isthe set of acquired 2D projection images, p is a parameter fixed so that∥I−ρ RR∥<1 and therefore guarantees convergence of the method and I isthe square unit matrix.

at each iteration, the algorithm V generates the sets of processedprojection images (p^((n)), q^((n)), t^((n)))=V[(p^((n-1)), q^((n-1)),t^((n-1)))] using the conjugate gradient method, typically given by:

$\quad\{ \begin{matrix}{\alpha^{({n - 1})} = {{- {\langle{q^{({n - 1})},{{RR}^{t}t^{({n - 1})}}}\rangle}}/{\langle{{{RR}^{t}t^{({n - 1})}},{{RR}^{t}t^{({n - 1})}}}\rangle}}} \\{p^{(n)} = {p^{({n - 1})} + {\alpha^{({n - 1})}t^{({n - 1})}}}} \\{q^{(n)} = {q^{({n - 1})} + {\alpha^{({n - 1})}{RR}^{t}t^{({n - 1})}}}} \\{\beta^{({n - 1})} = {{\langle{q^{(n)},{{RR}^{t}q^{(n)}}}\rangle}/{\langle{q^{({n - 1})},{{RR}^{t}q^{({n - 1})}}}\rangle}}} \\{{t^{(n)} = {q^{(n)} + {\beta^{({n - 1})}t^{({n - 1})}}}},}\end{matrix} $

in which p, q, t are such that p⁽⁰⁾ is arbitrary, typically zero orequal to s the vector of the acquired 2D projection images, andq⁽⁰⁾=t⁽⁰⁾=RR′p⁽⁰⁾−s.

at each iteration, the algorithm V_(θ(n)) generates a set of projectionimagesp^((n))=V_(θ(n))[p^((n-1))]=V_(θ(n))V_(θ(n-1))[p^((n-2))]=V_(θ(n))V_(θ(n-1)). . . V_(θ(1))[p⁽⁰⁾] with V_(θ(n))[p]=(I_(θ(n))−ρ_(θ(n))I_(θ(n))RR)p+ρ_(θ(n))I_(θ(n))s+I_(φ(n))p in which θ(n) and φ(n) are respectivelya First and a second sub-set of indices of acquired projectivemeasurements, defined for each iteration n, the two sub-sets beingseparate, their joining indexing the totality of the vector s of theacquired 2D projection images, and I_(θ(n)) and I_(φ(n)) are matriceswhose coefficients are all zero except for the indices of the diagonalbelonging to θ(n) and φ(n) respectively in which they equal 1, so thatp=I_(θ(n))+p+I_(θ(n))p, in which p⁽⁰⁾=ρ_(θ(0))I_(θ(0))s or else p⁽⁰⁾=0,and ρ_(θ(n)) is a fixe parameter so that ∥I_(θ(n))−ρ_(θ(n)) RR∥<1 andtherefore guarantees convergence of the method;

the method comprises extraction of a slice derived from thereconstructed 3D image;

the method comprises arbitrarily Fine adjustment of sampling of the 3Dimage.

According to a second aspect, the invention concerns a medical imagingsystem comprising means to implement a method according to the firstaspect of the invention.

According to a third aspect, the invention concerns a computer programmeproduct comprising programme code instructions to carry out the steps ofthe method according to the first aspect of the invention.

BRIEF DESCRIPTION OF THE DRAWINGS

Other characteristics and advantages of the invention will becomeapparent from the following description which is solely illustrative andnon-limiting and is to be read with reference to the appended drawingsin which, in addition to FIG. 1 already described:

FIG. 2 schematically illustrates a medical imaging system to implementthe method of the invention;

FIG. 3 schematically illustrates the steps of the method according tothe invention;

FIG. 4 schematically illustrates in detail a step of the method of theinvention.

DETAILED DESCRIPTION OF THE INVENTION Medical Imaging System

FIG. 2 schematically illustrates a medical imaging system 1 for theacquisition of 2D projection images to reconstruct a 3D image of anorgan.

Said system may be mammography apparatus for the detection andcharacterization of lesions in the screening, diagnosis and treatment ofbreast cancer.

The medical imaging system 1 comprises an image acquisition system 3, animage processing system 5 and a display system 4.

The acquisition system 3 is used to acquire a plurality of 2Dprojections of a region of interest—an organ—in a patient. Theacquisition system 3 notably consists of a detector D_(et) positionedfacing an X-ray source. The detector is a digital camera for example.The acquisition system is an acquisition system by X-rays, the systemcomprising any known means permitting the emission of X-rays onto theobject 2 and acquisition of the resulting images.

The display system 4 can be integrated in the image acquisition system 3or the image processing system 5, or it can be separated from theacquisition system 3 and the processing system 5. The display system 4is for example a computer screen, a monitor, flat screen, plasma screenor any type of known commercially available display device. The displaysystem 4 provides the practitioner with control over the reconstructionand/or display of the acquired 2D images.

The processing system 5 is adapted for implementation of processingmethods (e.g. reconstruction of a 3D image from 2D images). Theprocessing system 5 can be integrated in the image acquisition system 3or be separate from the image acquisition system 3. The processingsystem 5 is for example one or more computers, one or more processors,one or more microcontrollers, one or more micro-computers, one or moreprogrammable logic controllers, one or more specific applicationintegrated circuits, other programmable circuits or other devices whichinclude a computer such as a work station. The processing system 5 iscoupled with memory means 6 which can be integrated in or separate fromthe processing system 5. These means may consist of a hard disk or SSD,or any other removable re-write storage means (USB keys, memory cardsetc.). These memory means may serve to store a 3D image of the region ofthe organ viewed as an acquired or processed 2D image. They may be aROM/RAM memory of the processing system 5, a USB key, a memory card, acentral server memory. The processing system 5 may comprise a readdevice (not shown) e.g. a disk reader or CD-ROM reader to read theinstructions of the processing method (which is described below) from aninstruction medium (not shown) such as a disk or CD-ROM or moregenerally any removable memory medium or even via a network connection.As a variant, the processing system 5 may comprise a wire or wirelessnetwork connection device (not shown). As a variant, the processingsystem 5 carries out the instructions of the processing method(described below) stored in micro-software (not shown).

General Description of the Image Processing Method

The image processing method comprises an acquisition step S0 to acquirea plurality of 2D projection images s_(i) i=1, . . . , M>1, M being thenumber of acquired 2D projection images of an object, and a processingstep S1 to process the acquired 2D projection images.

It is the optimization of the processing step S1 which makes it possibleto limit the memory needs for storage of the 3D reconstruction.

More precisely, an iterative algorithm is applied to the acquiredprojection images to obtain processed 2D projection images p at eachiteration so that, at each iteration n, the 3D image of the object is alinear function of the processed 2D projection images.

In particular, the algorithm is defined by its iteration V whichgenerates a set of processed projection imagesp^((n))V[p^((n-1))]=V^(n)[p⁽⁰⁾] so that at each iteration n=1, . . . , Nthe 3D image of the object is a linear function of the processed 2Dprojection images, according to the property: f^((n))= Rp^((n)). It isthis equation which shows that, at each iteration n, it is possiblesimply to use a set of 2D projections to arrive at a reconstructed 3Dimage of the object.

It is noted that the processing step S1 comprises an initialization stepS10 and a step S11 during which iteration V is applied (see FIG. 4).

The processing step only requires storage of the sets of images(acquired and/or processed) instead of volumes.

Three algorithms are detailed below that are used in the 2D imagedomain. For each of these algorithms, their equivalence with a knowniterative reconstruction algorithm is shown used in the 3D image domain.

First Algorithm

A first conventional algorithm called <<successive approximations>>,consists of iteration A defined as follows:

A[f]=(I−ρ RR)f+ρ Rs.

in which I is the square unit matrix whose coefficients are all zeroexcept on the diagonal where they equal 1, R is a matrix such that RRfis a 3D image, and ρ an a priori fixed parameter so that ∥I−ρ RR∥<1 toguarantee convergence of the method.

It is known that the algorithm converges towards a solution of thesystem RRf= Rs. If R is reversible, this solution also meets Rf=s.

If R=R′, the estimate sequence converges towards a solution such that:

${\lim\limits_{narrow{(1)}}{A^{n}\lbrack f^{(0)} \rbrack}} = {f^{*} = {\min\limits_{f}{{{{Rf} - s}}^{2}.}}}$

In practice, a sufficiently large integer N is fixed for whichf^((n))=A[f^((n-1))] is computed for n=1, . . . , N with f⁽⁰⁾= Rp⁽⁰⁾with p⁽⁰⁾=0 or p⁽⁰⁾=ρs.

It is ascertained that in this standard algorithm, at each iteration, 3Dimages are handled (f, ρ RRf, ρ Rs, A[f]) which necessitates the storingof a high number of data (the entirety of at least one 3D image).

So as only to have 2D images to store, the algorithm defined byiteration V is applied, b computing the estimate sequencep^((n))=V[p^((n-1))] for n=1, . . . , N with p⁽⁰⁾=0 or p⁽⁰⁾=ρs, i.e.such that f⁽⁰⁾= Rp⁽⁰⁾.

Iteration V is defined for any vector p in space of the 2D projectionimages in the following manner:

V[p]=(I−ρR R )p+ρs.

It is noted that if f= Rp we have:

$\begin{matrix}{{\overset{\_}{R}( {V\lbrack p\rbrack} )} = {{{\overset{\_}{R}( {I - {\rho \; R\overset{\_}{R}}} )}p} + {\rho \overset{\_}{R}s}}} \\{= {{( {I - {\rho \; \overset{\_}{R}R}} )f} + {\rho \; \overset{\_}{R}s}}} \\{= {{A\lbrack f\rbrack}.}}\end{matrix}$

Therefore, with simple recurrence it has been shown that if f⁽⁰⁾= Rp⁽⁰⁾the two algorithms, one defined in the 2D image domain by iteration Vand the other in the 3D image domain by iteration A, give the samesuccession of 3D images f^((n)) since it has been shown that:

f ^((n)) =A ^(n) [f ⁽⁰⁾ ]= Rp ^((n)) = RV ^(n) [p ⁽⁰⁾].

It is to be pointed out, however, that this is not the case if there isthe constraint p=Rf since this would give:

$\begin{matrix}{{R( {A\lbrack f\rbrack} )} = {{{R( {I - {\rho \; R\overset{\_}{R}}} )}f} + {\rho \; R\overset{\_}{R}s}}} \\{= {{( {I - {\rho \; {\overset{\_}{R}}^{t}R}} )p} + {\rho \; R\; \overset{\_}{R}s}}} \\{\neq {{V\lbrack p\rbrack}.}}\end{matrix}$

Second Algorithm

In the particular case in which R=R′, the first algorithm—describedabove—gives us a solution to a quadratic optimization problem.

However, it is known that the algorithm of the conjugate gradient givesus a more efficient solution to this optimization problem, in that alargely reduced number of iterations needs to be computed to obtain avisually acceptable image.

On the other hand, from a storage viewpoint, the conjugate gradientrequires much more storage capacity than the first described algorithm.

It will be shown below that it can be implemented both in the 2Dprojection image domain and in the 3D image domain.

In manner known per se, in the 3D image domain, the algorithm of theconjugate gradient can be implemented by defining iteration A asfollows:

(f ^((n)) ,r ^((n)) ,d ^((n)))=A[(f ^((n-1)) ,r ^((n-1)) ,d ^((n-1)))]

with:

$\quad\{ \begin{matrix}{\alpha^{({n - 1})} = {{- {\langle{r^{({n - 1})},d^{({n - 1})}}\rangle}}/{\langle{{R^{t}{Rd}^{({n - 1})}},d^{({n - 1})}}\rangle}}} \\{f^{(n)} = {f^{({n - 1})} + {\alpha^{({n - 1})}d^{({n - 1})}}}} \\{r^{(n)} = {r^{({n - 1})} + {\alpha^{({n - 1})}R^{t}{Rd}^{({n - 1})}}}} \\{\beta^{({n - 1})} = {{\langle{r^{(n)},r^{(n)}}\rangle}/{\langle{r^{({n - 1})},r^{({n - 1})}}\rangle}}} \\{{d^{(n)} = {r^{(n)} + {\beta^{({n - 1})}d^{({n - 1})}}}},}\end{matrix} $

in which R is the matrix modelling the acquisition system and R′ thetransposed matrix, f^((n)) is an estimated 3D image of the solution tothe problem at iteration n, r^((n)) and d^((n)) are auxiliary 3D images,and <,> symbolizes the scalar product of two vectors so that α and β aretwo real numbers.

The three vectors in space of the 3D images are initialized as follows:

f ^((n)) =R′p ⁽⁰⁾ ,r ⁽⁰⁾ =d ⁽⁰⁾ =R′(Rf ⁽⁰⁾ −s)

in which s is the set of acquired 2D projection images, p⁽⁰⁾ isarbitrary but in practice is chosen to be equal to 0 or s.

It is simply ascertained that said implementation requires the storingof four volumes at the same time: f^((n)), r^((n)), d^((n)), R′Rd^((n)).

As previously, it is possible to dispense with the storage of these 3Dimages by changing the following variables:

$\begin{matrix}{{the}\mspace{14mu} {settings}\mspace{14mu} {for}} \\{{this}\mspace{14mu} {purpose}\mspace{14mu} {are}}\end{matrix}\{ \begin{matrix}p^{(n)} & {{such}\mspace{14mu} {that}} & {f^{(n)} = {R^{t}p^{(n)}}} \\q^{(n)} & {{such}\mspace{14mu} {that}} & {r^{(n)} = {R^{t}q^{(n)}}} \\t^{(n)} & {{such}\mspace{14mu} {that}} & {d^{(n)} = {R^{t}{t^{(n)}.}}}\end{matrix} $

in which (p^((n)), q^((n)), t^((n))) is a trio of vectors in space ofthe 2D projection images. The vector p⁽⁰⁾ is initialized as foriteration A and

q ⁽⁰⁾ =t ⁽⁰⁾ =RR′p ⁽⁰⁾ −s

Iteration V is defined as follows:

(p ^((n)) ,q ^((n)) ,t ^((n)))=V[(p ^((n-1)) ,q ^((n-1)) ,t ^((n-1)))]with

$\quad\{ \begin{matrix}{\alpha^{({n - 1})} = {{- {\langle{q^{({n - 1})},{{RR}^{t}t^{({n - 1})}}}\rangle}}/{\langle{{{RR}^{t}t^{({n - 1})}},{{RR}^{t}t^{({n - 1})}}}\rangle}}} \\{p^{(n)} = {p^{({n - 1})} + {\alpha^{({n - 1})}t^{({n - 1})}}}} \\{q^{(n)} = {q^{({n - 1})} + {\alpha^{({n - 1})}{RR}^{t}t^{({n - 1})}}}} \\{\beta^{({n - 1})} = {{\langle{q^{(n)},{{RR}^{t}q^{(n)}}}\rangle}/{\langle{q^{({n - 1})},{{RR}^{t}q^{({n - 1})}}}\rangle}}} \\{{t^{(n)} = {q^{(n)} + {\beta^{({n - 1})}t^{({n - 1})}}}},}\end{matrix} $

then, f^((N))=R′p^((N)) is determined, in which f is the 3D image of theobject to be reconstructed.

With this algorithm, it is indeed ascertained that at each iterationonly 2D images are stored, which makes it possible to obtain areconstruction that is less costly in terms of memory needs.

Therefore, it is possible to replace the intermediate steps f^((n)),r^((n)), d^((n)), R′Rd^((n)) in the 3D image domain by intermediatesteps in the 2D image domain: s^((n)), r^((n)), q^((n)), RR′q^((n)). Theequivalence is shown below between iteration A and iteration V of theconjugate gradient algorithm.

It is first shown that for any n=1, . . . , N, the scalar products α andβ are the same for iteration A in the 3D image domain as for iteration Vin the 2D image domain, since:

$\begin{matrix}{\alpha^{(n)} = {{- {\langle{r^{(n)},d^{(n)}}\rangle}}/{\langle{{R^{t}{Rd}^{(n)}},d^{(n)}}\rangle}}} \\{= {{\langle{{R^{t}q^{(n)}},{R^{t}t^{(n)}}}\rangle}/{\langle{{R^{t}{RR}^{t}t^{(n)}},{R^{t}t^{(n)}}}\rangle}}} \\{= {{- {\langle{q^{(n)},{{RR}^{t}t^{(n)}}}\rangle}}/{\langle{{{RR}^{t}t^{(n)}},{{RR}^{t}t^{(n)}}}\rangle}}}\end{matrix}$ and $\begin{matrix}{\beta^{(n)} = {{\langle{r^{({n + 1})},r^{({n + 1})}}\rangle}/{\langle{r^{(n)},r^{(n)}}\rangle}}} \\{= {{\langle{{R^{t}q^{({n - 1})}},{R^{t}q^{({n - 1})}}}\rangle}/{\langle{{R^{t}q^{(n)}},{R^{t}q^{(n)}}}\rangle}}} \\{= {{\langle{q^{({n + 1})},{{RR}^{t}q^{({n + 1})}}}\rangle}/{\langle{q^{(n)},{{RR}^{t}q^{(n)}}}\rangle}}}\end{matrix}$

Given that:

f^((n))=R′p^((n)),r^((n))=R′q^((n)),d^((n))=R′t^((n)),RR′t^((n)),RR′q^((n))

this leads to obtaining the following equivalences between iteration Vin the 2D image domain and iteration A in the 3D image domain:

$\quad\{ \begin{matrix}{{f^{(n)} = {R^{t}p^{(n)}}},{r^{(n)} = {R^{t}q^{(n)}}},{d^{(n)} = {R^{t}t^{(n\;)}}},{{RR}^{t}t^{(n)}},{{RR}^{t}q^{(n)}}} \\{q^{(n)}, {{RR}^{t}t^{(n)}}\Rightarrow\alpha^{(n)} } \\{p^{({n + 1})} = { {p^{(n)} + {\alpha^{(n)}t^{(n)}}}\Rightarrow f^{({n + 1})}  = {f^{(n)} + {\alpha^{(n)}d^{(n)}}}}} \\{q^{({n + 1})} = { {q^{(n)} + {\alpha^{(n)}{RR}^{t}t^{(n)}}}\Rightarrow r^{({n + 1})}  = {r^{(n)} + {\alpha^{(n)}R^{t}{Rd}^{(n)}}}}} \\ {{RR}^{t}q^{({n + 1})}}\Rightarrow\beta^{(n)}  \\{t^{({n + 1})} = { {q^{({n + 1})} + {\beta^{(n)}t^{(n)}}}\Rightarrow d^{({n + 1})}  = {r^{({n + 1})} + {\beta^{(n)}d^{(n)}}}}} \\{{{{RR}^{t}t^{({n + 1})}} =  {{{RR}^{t}q^{({n + 1})}} + {\beta^{(n)}{RR}^{t}t^{(n)}}}\Rightarrow\alpha^{({n + 1})} },}\end{matrix} $

we therefore obtain:

$\quad\{ \begin{matrix}{ p^{(0)}\Rightarrow f^{(0)}  = {R^{t}p^{(0)}}} \\{q^{(0)} = { {{t^{(0)}{RR}^{t}p^{(0)}} - s}\Rightarrow r^{(0)}  = {d^{(0)} = {R^{t}( {{Rf}^{(0)} - s} )}}}} \\{ {{determination}\mspace{14mu} {of}\mspace{14mu} {RR}^{t}t^{(0)}}\Rightarrow\alpha^{(0)} ,{{{RR}^{t}q^{(0)}} = {{RR}^{t}t^{(0)}}}} \\{p^{(1)} = { {p^{(0)} + {\alpha^{(0)}t^{(0)}}}\Rightarrow f^{(1)}  = {f^{(0)} + {\alpha^{(0)}d^{(0)}}}}} \\{{q^{(1)} = { {q^{(0)} + {\alpha^{(0)}{RR}^{t}t^{(0)}}}\Rightarrow r^{(1)}  = {r^{(0)} + {\alpha^{(0)}R^{t}{Rd}^{(0)}}}}} {{determination}\mspace{14mu} {of}\mspace{14mu} {RR}^{t}q^{(1)}}\Rightarrow\beta^{(0)} } \\{t^{(1)} = { {q^{(1)} + {\beta^{(0)}t^{(0)}}}\Rightarrow d^{(1)}  = {r^{(1)} + {\beta^{(0)}d^{(0)}}}}} \\{{{RR}^{t}t^{(1)}} =  {{{RR}^{t}q^{(1)}} + {\beta^{(1)}{RR}^{t}t^{(1)}}}\Rightarrow{\alpha^{(1)}.} }\end{matrix} $

We have shown, by recurrence, the equivalence between iteration V in the2D image domain and iteration A in the 3D image domain.

Finally, as for the preceding algorithm, the 3D image of the object isobtained by applying the linear constraint, here f^((N))=R′p^((N)).

The formula of the conjugate gradient given above is that of Fletcherand Reeves, but similar reasoning is possible with other formulations ofthe conjugate gradient (Polak-Ribiere, Graham-Schmidt, . . . )

Similarly, the result applies for any matrix R or if the product RR is apositively defined symmetrical matrix.

Third Algorithm

This third algorithm called “block iterative” generalizes the first inthat all the indices of the measurements will be split into two sub-setsof indices.

The first sub-set is denoted θ and the second is denoted (1). The twosub-sets are separate, and their joining indexes all the 2D projectionimages.

In a first variant, the sub-set θ contains a single index (i.e. a singlemeasurement point of the detector for a single given acquisition angle)whereas φ indexes all the other measurements (i.e. all except one). Thispartitioning is used in the algorithm “Algebraic ReconstructionTechnique (ART)”.

In a second variant, the sub-set θ indexes all the 2D images for asub-set of angles at which the projection images are acquired, thesub-set φ indexing the complementary acquired angle positions. Thispartitioning is used in the algorithm “Simultaneous AlgebraicReconstruction Technique (SART)”. It is noted that the case in which θcontains all the indices and φ is empty corresponds to the case of thefirst algorithm.

I_(θ) (respectively I_(φ)) is used to denote the matrix whosecoefficients are all zero except for the indices of the diagonalbelonging to θ (respectively φ) in which they equal 1. Since (θ, φ) is apartition, the unit matrix in space of the 2D images is I_(θ)+I_(φ).

It is noted that the product I_(θ)R (respectively I_(φ)R) is therestriction of R at the indices of θ (respectively φ). Similarly, it isnoted that the product RI_(θ) (respectively RI_(φ)) is the restrictionof R at the indices of θ (respectively φ).

On the basis of these definitions, the algorithm in the 3D image domaincan be implemented by iteration 4, defined as follows:

A _(θ) [f]=(I−ρ _(θ RI) _(θ) R)f+ρ _(θ) RI _(θ) s

in which 1 is the unit matrix in the 3D image domain and ρ_(θ) an apriori fixed parameter so that μI−ρ_(θ) RI_(θ)R∥<1.

Evidently, any reconstruction uses all the measurements to determine asolution, therefore the partitioning (θ, φ) of indices is not constantbut changes with each iteration. It is therefore a plurality ofiterations A_(θ).

The algorithm then generates an estimate sequence f^((n)) such that:

f ^((n)) =A _(θ(n)) [f ^((n-1)) ]=A _(θ(n)) A _(θ(n-1)) [f ^((n-2) ]=A_(θ(n)) A _(θ(n-1)) . . . A _(θ(1)) [f ⁽⁰⁾]

It is noted that at each iteration n, the estimate is obtained from arestricted set of measurements θ(n), so that the computing time of theiteration is reduced. It is known that these techniques thereforeprovide solutions to the problem raised with reduced computing time. Onthe other hand, the problem associated with the storing of 3D imagesremains unchanged and independent of the chosen partitioning (θ, φ).

It is shown here in novel manner that it is possible as previously todetermine iterations V_(θ) in the 2D projection image domain which willlead to the same estimates as those obtained with iterations A_(θ).

So as only to have 2D images to store, the algorithm defined byiteration V_(θ) is applied by computing the estimate sequencep^((n))=V_(θ(n))[p^((n-1))] for n=1, . . . , N with p⁽⁰⁾ such thatf⁽⁰⁾=R′p⁽⁰⁾ and in accordance with the same θ(n) partitioning as used inequivalent iteration A_(θ). Iteration V_(θ) is defined for any vector pin space of the 2D projection images and for any (θ, φ) partitioning ofmeasurements such that p=I_(θ)p+I_(φ)p, in the following manner:

V _(θ) [p]=(I _(θ)−ρ_(θ) I _(θ) R R )p+ρ _(θ) I _(θ) s+I _(φ) p.

Finally, as for the two other algorithms, the 3D image of the object isobtained by determining f^((N))= Rp^((N)).

This third algorithm—as can be ascertained—at iteration n only modifiesthe variables indexed by θ.

It is shown that, irrespective of (θ, φ) partitioning, if f= Rp, it ispossible to change over from the algorithm defined by iterations A_(θ)in the 3D image domain, to the algorithm defined by iterations V_(θ) inthe 2D image domain, in the following manner:

$\begin{matrix}{{\overset{\_}{R}( {V_{\Theta}\lbrack p\rbrack} )} = {{\overset{\_}{R}( {{( {I_{\Theta} - {\rho_{\Theta}I_{\Theta}\; R\overset{\_}{R}}} )p} + {\rho_{\Theta}I_{\Theta}s}} )} + {\overset{\_}{R}I_{\Theta}p}}} \\{= {{\overset{\_}{R}( {{I_{\Theta}p} + {I_{\Theta}p}} )} - {\rho_{\Theta}\overset{\_}{R}I_{\Theta}R\overset{\_}{R}p} + {\rho_{\Theta}\; \overset{\_}{R}I_{\Theta}s}}} \\{= {f - {\rho_{\Theta}\overset{\_}{R}I_{\Theta}{Rf}} + {\rho_{\Theta}\overset{\_}{R}I_{\Theta}s}}} \\{= {A_{\Theta}\lbrack f\rbrack}}\end{matrix}$

As previously, iteration V_(θ), at any time (initialization n=0,intermediate computing 0<n<N archiving n=N) avoids having to store theentirety of a 3D image whose size is much larger than all the acquiredmeasurements.

Application to Few-View Tomography

As already mentioned, one problem in tomography is the low number ofprojections when it is not possible to rotate completely around theobject it is desired to image, or if the acquisition must be sub-sampledto limit the examination time or the X-ray dose to the patient. As aresult, the size of the 3D image it is desired to reconstruct is muchgreater than the size of all the acquired 2D projection images.

By using the reconstruction algorithms such as presented, no volume isstored which means that the sampling of the volume can be finelyadjusted and arbitrarily without any impact on memory storage needs.

Since the final step f^((N))= Rp^((N)) which, for only one slice to beviewed, is merely a linear operator applied to a low number of processedacquired 2D projection images, it can be performed in real time whenviewing the slice without ever having to store all the slices of thevolume.

It is noted that, depending on the desired resolution for the 3D imageto be reconstructed, the number of iterations can be increased withoutincreasing memory storage needs.

The methods presented are particularly suitable for GPUs (GraphicProcessor Units) which have a good performance level with respect tocomputing time, but which have limited Memory.

1. A method for processing images obtained by tomography ortomosynthesis, the method comprising: acquiring a plurality of 2Dprojection images of an object using an imaging system, the acquiring aplurality of 2D projection images of an object being defined by Rf=swherein s is a vector of the acquired projections, R is a projectionoperator which models the imaging system and f is the 3D image of theobject to be reconstructed with knowledge of R and s; and processing theacquired 2D projection images, wherein the processing the acquired 2Dimages comprises applying an iterative algorithm defined by itsiteration, wherein at least one set of processed projection images isgenerated at each iteration, the iteration being defined so that at eachiteration the 3D image of the object is a linear function of theprocessed 2D projection images according to the property f^((n))=Rp^((n)) wherein R is a matrix such that RRf is a 3D image.
 2. Themethod according to claim further comprising reconstructing the 3D imageof the object to be determined according to the property f^((N))=Rp^((N)), wherein N is the number of fixed iterations.
 3. The methodaccording to claim 1, wherein at each iteration, the algorithm generatesa set of processed projection images, the iteration being defined forany vector in space of the 2D projection images by:V[p]=(I−ρR R )P+ρs, wherein p⁽⁰⁾=ρs or p⁽⁰⁾=0, s is the set of acquired2D projection images, ρ is a fixed parameter so that ∥I−ρ RR∥<1 andthereby guarantees convergence of the method, and I is the square unitmatrix.
 4. The method according to claim 1, wherein at each iterationthe algorithm generates sets of processed projection images as per theconjugate gradient method: $\quad\{ \begin{matrix}{\alpha^{({n - 1})} = {{- {\langle{q^{({n - 1})},{{RR}^{t}t^{({n - 1})}}}\rangle}}/{\langle{{{RR}^{t}t^{({n - 1})}},{{RR}^{t}t^{({n - 1})}}}\rangle}}} \\{p^{(n)} = {p^{({n - 1})} + {\alpha^{({n - 1})}t^{({n - 1})}}}} \\{q^{(n)} = {q^{({n - 1})} + {\alpha^{({n - 1})}{RR}^{t}t^{({n - 1})}}}} \\{\beta^{({n - 1})} = {{\langle{q^{(n)},{{RR}^{t}q^{(n)}}}\rangle}/{\langle{q^{({n - 1})},{{RR}^{t}q^{({n - 1})}}}\rangle}}} \\{{t^{(n)} = {q^{(n)} + {\beta^{({n - 1})}t^{({n - 1})}}}},}\end{matrix} $ wherein p, q, t are such that p⁽⁰⁾ is arbitrary,typically zero or equal to the vector of the acquired 2D projectionimages, and en q⁽⁰⁾=t⁽⁰⁾=RR′p⁽⁰⁾−s.
 5. The method according to claim 1,wherein at each iteration the algorithm generates a set of projectionimages withV _(θ) [p]=(I _(θ)−ρ_(θ) I _(θ) R R )p+ρ _(θ) I _(θ) s+I _(φ) p, whereinθ(n) and φ(n) are respectively a first and a second sub-set of indicesof acquired projective measurements, defined for each iteration, the twosub-sets being separate, their joining indexing the entirety of thevector of the acquired 2D projection images, and I_(θ(n)) and I_(φ(n))are matrices whose coefficients are all zero except for the indices ofthe diagonal belonging to θ(n) and φ(n) respectively in which theirvalue is 1, so that p=I_(θ(n))p+I_(φ(n))p, in whichp⁽⁰⁾=ρ_(θ(0))I_(θ(0))s or else p⁽⁰⁾=0, and ρ_(θ(n)) is a fixed parameterso that ∥I_(θ(n))−ρ_(θ(n)) RR∥<1 and thereby guarantees the convergenceof the method.
 6. The method according to claim 1, further comprisingextracting a slice derived from the reconstructed 3D image.
 7. Themethod according to claim 1, further comprising arbitrarily adjustingthe sampling of the image.
 8. A medical imaging system comprising: anacquisition system for the acquisition of a plurality of 2D projectionimages of an object, comprising a radiation source and a sensor; astorage means configured to store acquired 2D projection images; and aprocessing system configured to: acquire a plurality of 2D projectionimages of an object using an imaging system, the acquiring a pluralityof 2D projection images of an object being defined by Rf=s, wherein s isa vector of the acquired projections, R is a projection operator whichmodels the imaging system and f is the 3D image of the object to bereconstructed with knowledge of R and s; and process the acquired 2Dprojection images, wherein the processing the acquired 2D imagescomprises applying an iterative algorithm defined by its iteration,wherein at each iteration n=1, . . . , N generates at least one set ofprocessed projection images, the iteration being defined so that at eachiteration the 3D image of the object is a linear function of theprocessed 2D projection images according to the property f^((n))=Rp^((n)), wherein R is a matrix such that RRf is a 3D image.